% master_simulation_E.m 
clear all;
clc;
close all;
tStart = tic;
logfilename = 'LOG_MASTER_SIMULATION.txt';
if (exist(logfilename)) 
    delete(logfilename);
end
diary(logfilename);
alphavec=[-0.1:0.005:0.3]; 
betavec=[-0.1:0.005:0.3];
n=size(alphavec,2);
alphax=0.103; 
betax=0.213; 
k=0.4;      
Kmat=[1 k;k 1]; 
% ***********************************
% **** 1. baseline equilibrium  *****
% ***********************************
rmat=zeros(n,n); 
lmat=zeros(n,n); 
hmat=zeros(n,n); 
BB=ones(2,1); 
AA=ones(2,1);
for i=1:n
    alpha=alphavec(i); 
    for j=1:n
        beta=betavec(j); 
        theta=0.9; R_i=[0.5;0.5]; L_i=[0.5;0.5]; Q_i=[1;1]; 
        m_equivsolve_level; 
        rmat(i,j)=R(1)/R(2); lmat(i,j)=L(1)/L(2); hmat(i,j)=H(1)/H(2); 
    end 
end 

m_popfigure; print(gcf, '-dpdf', '../../output/appendixE/population_dynamics_base.pdf','-fillpage');
m_empfigure; print(gcf, '-dpdf', '../../output/appendixE/employment_dynamics_base.pdf','-fillpage');

% ********************************************
% ***** 2. different equilibrium *************
% ********************************************
rmat=zeros(n,n); 
lmat=zeros(n,n); 
hmat=zeros(n,n); 
BB=ones(2,1); 
AA=ones(2,1);
for i=1:n
    alpha=alphavec(i); 
    for j=1:n
        beta=betavec(j); 
        theta=0.9; R_i=[0.8; 0.2]; L_i=[0.8; 0.2]; Q_i=[1;1]; 
        m_equivsolve_level; 
        rmat(i,j)=R(1)/R(2); lmat(i,j)=L(1)/L(2); hmat(i,j)=H(1)/H(2); 
    end 
end 

m_popfigure_2; print(gcf, '-dpdf', '../../output/appendixE/population_dynamics_multiple.pdf','-fillpage');
m_empfigure_2; print(gcf, '-dpdf', '../../output/appendixE/employment_dynamics_multiple.pdf','-fillpage');

%% **************************************************
% ******** 3. different parameter (theta) **********
% **************************************************
rmat=zeros(n,n); 
lmat=zeros(n,n); 
hmat=zeros(n,n); 
BB=ones(2,1); 
AA=ones(2,1);
for i=1:n
    alpha=alphavec(i); 
    for j=1:n
        beta=betavec(j); 
        theta=0.5; R_i=[0.8;0.2]; L_i=[0.8;0.2]; Q_i=[1;1];  
        m_equivsolve_level; 
        rmat(i,j)=R(1)/R(2); lmat(i,j)=L(1)/L(2); hmat(i,j)=H(1)/H(2); 
    end 
end 

m_popfigure; print(gcf, '-dpdf', '../../output/appendixE/population_dynamics_lowtheta.pdf','-fillpage');
m_empfigure; print(gcf, '-dpdf', '../../output/appendixE/employment_dynamics_lowtheta.pdf','-fillpage');

%% **********************************************************************************
% ******** 4. different fundamentals  (location 2 has small advantage) **************
% ***********************************************************************************
rmat=zeros(n,n); 
lmat=zeros(n,n);
hmat=zeros(n,n); 
BB=[1.05;1]; 
AA=[1; 1];
for i=1:n
    alpha=alphavec(i); 
    for j=1:n
        beta=betavec(j); 
        theta=0.9; R_i=[0.5;0.5]; L_i=[0.5;0.5]; Q_i=[1;1]; 
        m_equivsolve_level; 
        rmat(i,j)=R(1)/R(2); lmat(i,j)=L(1)/L(2); hmat(i,j)=H(1)/H(2); 
    end 
end 

m_popfigure_2; print(gcf, '-dpdf', '../../output/appendixE/population_dynamics_fundadvantage_norecovery.pdf','-fillpage');
m_empfigure_2; print(gcf, '-dpdf', '../../output/appendixE/employment_dynamics_fundadvantage_norecovery.pdf','-fillpage');

rmat=zeros(n,n); lmat=zeros(n,n); hmat=zeros(n,n); BB=[1.05; 1]; AA=[1;1];
for i=1:n
    alpha=alphavec(i); 
    for j=1:n
        beta=betavec(j); 
        theta=0.9; R_i=[0.8;0.2]; L_i=[0.8;0.2]; Q_i=[1;1]; 
        m_equivsolve_level; 
        rmat(i,j)=R(1)/R(2); lmat(i,j)=L(1)/L(2); hmat(i,j)=H(1)/H(2); 
    end 
end 
 
m_popfigure_2; print(gcf, '-dpdf', '../../output/appendixE/population_dynamics_fundadvantage_recovery.pdf','-fillpage');
m_empfigure_2; print(gcf, '-dpdf', '../../output/appendixE/employment_dynamics_fundadvantage_recovery.pdf','-fillpage');
disp(['Elapsed time is: ', num2str(toc / 60), ' minutes']);
diary off; 